Comparison of ChangeDetection Techniques for Monitoring Tropical Forest Clearing and Vegetation Regrowth in a Time Series

The once remote and inaccessible forests of Guatemala's Maya Biosphere Reserve (MBR) have recently experienced high mtes of deforestation corresponding to human migration and expansion of the agricultural frontier. Given the importance of land-cover and land-use change data in conservation planning, accurate and efficient techniques to detect forest change from multi-tempoml satellite imagery were desired for implementation by local conservation organizations. Three dates of Landsat Thematic Mapper imagery, each acquired two years apart, were radiometrically normalized and preprocessed to remove clouds, water, and wetlands, prior to employing the change-detection algorithm. Three changedetection methods were evaluated: normalized difference vegetation index (NDVI) image differencing, principal component analysis, and RGB-NDVI change detection. A technique to generate reference points by visual interpretation of color composite Landsat images, for Kappa-optimizing thresholding and accuracy assessment, was employed. The highest overall accuracy was achieved with the RGB-IVDVI method (85 percent). This method was also preferred for its simplicity in design and ease in interpretation, which were important considerations for transferring remote sensing technology to local and international non-governmental organizations. Introduction With rapid changes in land-cover occurring over large areas, remote sensing technology is an essential tool in monitoring tropical forest conditions. The remote and inaccessible nature of many tropical forest regions limits the feasibility of groundbased inventory and monitoring methods for extensive land areas. Initiatives to monitor land-cover and land-use change are increasingly reliant on information derived from remotely sensed data. Such information provides the data link to other techniques designed to understand the human processes behind deforestation (Lambin, 1994; Rindfuss and Stern, 1998). An array of techniques are available to detect land-cover changes from multi-temporal remote sensing data sets Uensen, 1996; Coppin and Bauer, 1996). The goal of change detection is to discern those areas on digital images that depict change features of interest (e.g., forest clearing or land-covert land-use change) between two or more image dates. One method, image differencing, is simply the subtraction of the pixel digital values of an image recorded at one date from the corresponding pixel values of the second date. The histogram Maint: Imago Allalysis Lnl~oratory. L)cparlment o f E'oresl Manag(:ment. 5755 Nutting Hall, Il~iive~.silv of Maine. Orono, ME 04469 [~l11ayr!s@imc~1~fa.111:1i11e.e(1~1; sadr:r@~~rnonf~~.~~i~~i~~e~r:~l~~). of the resulting image depicts a range of pixel values from negative to positive numbers, where those clustered around zero represent no change and those at either tail represent reflectance changes from one image date to the next (Jensen, 1996). This method has been documented widely in changedetection research (Singh, 1986; Muchoney and Haack, 1994; Green et al., 1994; Coppin and Bauer, 1996; Macleod and Congalton, 1998). Some investigators favor this method for its accuracy, simplicity in computation, and ease in interpretation. One difficulty encountered in employing image differencing for change detection is the selection of the appropriate threshold values in the histogram that separates real and spurious change. The subjectivity of threshold placement may be improved by the analyst's familiarity with the study area as well as access to ancillary data such as field information, GIS data, and/or matching dates of aerial photography. Fung and LeDrew (1988) tested quantitative methods for developing these threshold levels using accuracy indices. They recommended the Kappa coefficient of agreement in determining an optimal threshold level, being based on an error matrix of image data against known reference data. Image differencing, although mathematically simple, allows for only one band of information to be processed at a time. Other techniques incorporate multiple bands of data for change detection. Several studies have demonstrated the utility of the principal component analysis (PCA) technique in multi-temporal image analysis (Byrne et al., 1980; Fung and LeDrew, 1988; Muchoney and Haack, 1994; Coppin and Bauer, 1996; Macleod and Congalton, 1998). The results of using the PCA transform on two dates of imagery are contrary to that of its typical, one-date transformations. In multi-temporal analysis, the first two components tend to represent variation associated with unchanged land-cover and overall image noise (i.e., atmospheric and seasonal variation), while the third and later components are of more interest in identifying change areas (Byrne et al., 1980). Previous studies have confirmed that the minor components have been successful in detecting landcover changes (Byrne et al., 1980; Fung and LeDrew, 1987) when the areas affected by change of interest occupy a small proportion of the study area (Fung and LeDrew, 1987; Macleod and Congalton, 1998). Image differencing using band ratios or vegetation indices is another technique commonly employed for land-cover I'hotogra~nm(!tric: Ellgil~cering & Ktr~l~ote Sensing Vol. 67. No. $1, Ssptc:nlt)er 2001. 111). 10(i7-1075. ~ 0 ~ ~ ~ 1 1 1 2 / 0 1 / 6 7 0 ~ ~ 1 0 t i 7 S 3 . 0 0 / 0 6? 2001 A~nc!ric:an Society for Pirotogr~rmmctrv H I I ~ lio1110t~: Sensing change detection. For example, the normalized difference vegetation index (NDVI) was developed for use in identifying health and vigor in vegetation, as well as for estimates of green biomass. The NDvI, the normalized difference of brightness values from the near infrared and visible red bands, has been found to be highly correlated with crown closure, leaf area index, and other vegetation parameters (Tucker, 1979; Sellers, 1985; Singh, 1986; Running et al., 1986). Lyon et al. (1998) compared seven vegetation indices to detect land-cover change in a Chiapas, Mexico study site. They reported that the NDVI was least affected by topographic factors and was the only index that showed histograms with normal distributions. Change in canopy cover or vegetation biomass can be detected by analyzing NDVI values from separate dates (e.g., NDvI image differencing). Sader and Winne (1992) developed a technique to visualize change using three dates of NDVI imagery concurrently and interpretation concepts of color additive theory. By simultaneously projecting each date of NDVI through the red, green, and blue (RGB) computer display write functions, major changes in NDVI (and, hence, green biomass) between dates will appear in combinations of the primary (RGB) or complimentary (yellow, magenta, cyan) colors. Knowing which date of NDVI is coupled with each display color, the analyst can visually interpret the magnitude and direction of biomass changes in the study area over the three dates. Automated classification can be performed on three or more dates of NDVI by unsupervised cluster analysis (Sader et al., 2001). Change and no-change categories are labeled and dated by interpreter analysis of the cluster statistical data and guided by visual interpretation of R G B ~ v I color composites. Study Area and Background Spanning approximately 2 million hectares of northern Guatemala, the Maya Biosphere Reserve (MBR) is an area of lowland tropical forests and expansive freshwater wetlands, part of the largest continuous tropical moist forest remaining in Central America (Nations et al., 1998). The MBR is a complex of delineated management units, including five national parks, four biological reserves (biotopos), a multiple use zone, and a buffer zone (Figure I). The once remote and inaccessible forests of the region have experienced high rates of deforestation in the last decade, corresponding to human migration and expansion of the agricultural frontier (Sader et al., 1997). Sader and colleagues (Sader et al., 1997; Sader et al., 2001) have monitored rates and trends of forest clearing using Landsat Thematic Mapper (TM) imagery from the mid-1980s to late 1990s. Guatemalan government agencies and non-governmental organizations (NGOS) rely on regularly updated maps of the MBR to monitor deforestation patterns and disturbance in sensitive areas of the reserve. International donor agencies require the NGOs to quantify forest clearing rates at two-year intervals. Accurate and efficient techniques for extracting quantitative forest-change data from remotely sensed images are needed to support the MBR forest monitoring program. Furthermore, these data are needed for analysis with community level socio-economic survey data concerning the driving forces of environmental change in the MBR (Schwartz, 1998; Hayes, 1999). This paper describes the techniques used to process and validate multi-temporal Landsat TM imagery (three dates) for obtaining time-series forest clearing and regrowth data in the MBR. Three change-detection methods are compared: NDM image differencing, PCA change detection, and RGB-NDVI classification. A visual interpretation technique to generate reference points from color composite Landsat images, for selecting Kappa-optimizing thresholds and for assessment of classification accuracy, is described. The goal is to determine the most accurate and efficient method to detect forest change in the 1068 September 2001 Figure 1. Location of the study area (Landsat WRS Path 20/Row 48, 1997 TM band 5 shown) in relation to the management units of the Maya Biosphere Reserve, El Peten, Guatemala. MBR'S tropical moist forest and to facilitate the transfer of this technology to the local NGOs. Data Acquisltlon and Pre-Processing Three dates of Landsat TM imagery (1993, 1995, 1997) for Worldwide Reference System path 20, row 48 were acquired. This Landsat scene comprises approximately 90 percent of the MBR and buffer zone (Figure 1). To reduce scene-to-scene variation due to sun angle, soil moisture, atmospheric condition, and vegetation phenology differences, all data were collected between the months of March and May, corresponding to the MBR'S dry season. Each scene was georeferenced to a previously rectified 1995 TM image. TM bands 3 (visible red), 4 (near inhared), and 5 (mid-infrared) were extracted from the original TM data sets to reduce between-band correlation, data volume, and processing time. Previous studies have shown that selecting one band each from the visible, near infrared, and mid-infrared spectral regions results in the optimal waveband combination for vegetation discrimination (DeGloria, 1984; Horler and Ahern, 1986; Sader, 1989). Bands 3,4 , and 5 were input into "isodata" (ERDAS, 1997), an unsupervised classification module, to produce 200 spectral clusters. Binary images were created to isolate water, clouds, and cloud shadows through a combination of analyst definition of cloudlwater clusters and on-screen editing. A previously developed image of non-forested wetlands and natural savannas was also added PHOTOGRAMMETRIC ENGINEERING & REMOTE SENSING to the cloud and water image. These classes, being of no interest to forest clearing and regrowth analysis, were masked for all and dates of imagery to avoid confusion in the change-detection classification. Radiometric Nonnallzatlon A relative radiometric calibration technique was applied to each band from each date of imagery. The technique incorporated linear regression methods reported by Eckhardt et al. (1990), Hall et al. (1991), and Jensen et al. (1995). The 1997 TM scene, which was corrected for sensor gain and bias, was used as the reference image to which the 1993 and 1995 data were normalized. First, normalization targets were selected from the wet (e.g., deep, clear water) and dry (e.g., urban features) nonvegetated extremes of each band (TM 3,4, and 5) at each date (1993,1995,1997) by visual interpretation of the imagery and querying the digital numbers of pixels representing these features. The selection criteria were based on procedures outlined by Eckhardt et al. (1990). Each target consisted of an analystdefined area of interest (AOI), which included the greatest number of pixels covering the target, whose digital numbers (DNS) were located at the extremes of the image histogram and collectively contained low variance. The mean value of the pixel DNs was generated for each of the normalization target AOIS (each band, each date). The parameters used in the linear regression equation were calculated by the following "rectification transform" (Hall et al., 1991): where Br is the mean DN for the bright target of the reference image, Bs is the mean DN for the bright target of the subject image, Dr is the mean DN for the dark target of the reference image, and Ds is the mean DN for the dark target of the subject image. Using linear regression, the corrected pixel values for the subject image (Y) were calculated from the original DN (XI, for each band (11, by the following equation: Changehtectlon Methods Three change-detection methods (NDVI differencing, PCA, and RGB-NDW classification) were independently applied to the cloudlwater-masked and radiometrically normalized timeseries Tbi data set. A three-date forest change-detection classification of the selected study area was generated from each method. Each method was evaluated and compared with the other methods on its ability to classify temporal states in forest cover (i.e., cleared, regrown, no change) over the three time periods. The methods were evaluated and contrasted on the basis of classification accuracy (Congalton, 1991), efficiency in computation and processing, and ease in interpretation. NDVI Image Differencing Difference images were created by first calculating NDvI values for each date ( j ) of imagery by the following equatioa: Principal Component Analysis The principal component transformation was performed separately on two data sets (1993 and 1995,1995 and 1997) using three TM bands (3,4, and 5) for each date. Each two-date data set contained six bands. The transformation used the "prince" routine (ERDAS 1997), modified to calculate the transform from a correlation matrix of the data set. Several authors have compared this "standardized" approach to PCA against transformations based on the covariance matrix (Conese et al., 1988; Eastman and Falk, 1993; Rencher 1995). Reported advantages of the standardized approach include improved interpretability, the isolation of seasonal effects and variability due to noise, better statistical control, and more precise classification. For each data set, the "standardized" PCA routine output included six component images, a table of eigenvalues quantifying the proportion of variance explained by each component, and a matrix of eigenvectors (weights or factor loadings) depicting between-date correlation for each band with each component. Components that represent change typically show an ibsence of correlation amone bands between dates IBvrne et al.. 19801. The component that'best highlights the chGge of interest is ' chosen for thresholding, using visual interpretation of component images and analysis of the eigenvector matrix. Image interpretation was based on the assessment of spatial continuity, by seeking out the components that express the differences in the changes of interest as spatially discontinuous areas within the image. The eigenvector analysis examined the algebraic signs on the weights. Differences between dates are expressed by the weight of one band at one date having an opposite sign to that of the same band of the other date. Based on these criteria, two of the six components (components 3 and 4 for each two-date data set) were selected from the PCA for thresholding of no-change and change areas. Of these two components, the one that showed the highest ability to threshold forest clearinglno-changelregrowth (i.e., the highest estimated Kappa according to the reference sample points) was chosen for final classification. RGB-NDVZ Classification NDVI values from three dates (as calculated by Equation 3) were classified into 50 spectral clusters. For each cluster class, the mean NDVI values at each date (1993,1995,1997) were categorized as very high, high, medium-high, medium, medium-low, low, or very low, based on the distribution of ~ ~ v r values over the study area. These levels of NDVI were established on the observation that, because most of the study area is composed of undisturbed forest, values within r 0.5 standard deviations from the mean represented high green biomass (high mean NDVI). The other NDVI levels were set at intervals of 0.5 standard deviations outward from the mean. Each cluster was examined for changes in NDVI levels over time. Clusters were named according to type of change (clearing, regrowth, or no change) and the corresponding time period(s) of change according to the NDW levels as they related to three-date RGBNDVI interpretation (Plate 1). (TM4 m 3 ) NDw[jl = (TMI + T M ~ ) (3) Classlfylng the Change Images Both the NDVI differencing and the PCA methods result in images with an 8-bit (0 to 255) data range. Thresholds must be Two difference images were created by subtracting one date of identified along the histograms to separate change (both clearNDVI values from those of the previous date, so that ing and regrowth) from no change. Threshold levels were set PHOTOGRAMMETRIC ENGINEERING 81 REMOTE SENSING September 200 / 1089